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Abstract. In this paper we numerically solve the eigenvalue problem Au+Xu — on the fractal re- 
gion defined by the Koch Snowflake, with zero-Dirichlet or zero-Neumann boundary conditions. The 
Laplacian with boundary conditions is approximated by a large symmetric matrix. The eigenvalues 
and eigenvectors of this matrix are computed by ARPACK. We impose the boundary conditions in 
a way that gives improved accuracy over the previous computations of Lapidus, Neuberger, Renka 
& Griffith. We extrapolate the results for grid spacing h to the limit h — > in order to estimate 
eigenvalues of the Laplacian and compare our results to those of Lapdus et al. We analyze the 
symmetry of the region to explain the multiplicity-two eigenvalues, and present a canonical choice 
' of the two eigenfunctions that span each two-dimensional eigenspace. 

6 

in 

1. Introduction. 

| In this paper we approximate solutions to the two eigenvalue problems 

^ ■ Au + Xu = in Q Au + Xu = in 

( L1 ) « = on3!](D) ^ = on 90 (N), 

or] 

where A is the Laplacian operator, and f! C M 2 is the (open) region whose boundary dQ is the Koch 
snowflake. For convenience, we refer to O as the Koch snowflake region. The boundary conditions 
are zero-Dirichlet, or zero-Neumann, respectively. 

These boundary value problems must be interpreted in the variational sense (see |Lapidus, 1991| ) 
C"-- but we avoid the subtleties of functional analysis by discretizing the problem. We use a triangular 

grid of points to approximate the snowflake region. Then, we identify u : O — >• M with u G M. N , 
where N is the number of grid points in Q. That is, 

u(xi) fa Ui 

at grid points Xi 6 M 2 , i £ {1,2,3, . . . ,N}. The discretized Laplacian is the symmetric matrix L, 
with the property 

N 

(-Au)(xi) fa (Lu)i = y~]LjjUj. 

j=i 

Of course, a specific grid and a scheme for enforcing the boundary conditions are needed to define 
L. This is described in Section 2. Then, the eigenvalues and eigenfunctions of L approximate the 
eigenvalues and eigenfunctions defined by (jl.ip . The eigenvalues and eigenvectors of L are our 
approximations of the eigenvalues < Ai < A2 < A3 < ••• < X^ ■ ■ ■ — > 00 and the corresponding 
eigenfunctions {V'fclfcLi of the negative Laplacian —A. 

The Koch snowflake is a well known fractal, with Hausdorff dimension log 3 4. Following Lapidus, 
Neuberger, Renka, and Griffith |Lapidus et al., 19 96 , we take our snowflake to be inscribed in a 



> 

in 



X 



circle of radius ^ centered about the origin. With this choice, the polygonal approximations used 
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in the fractal construction have side length that are powers of 1/3. In |Lapidus et al., 1996| , a 
triangular grid with spacing h = /ilnr(^) = 1/3 was used to approximate the eigenfunctions. Here 
£ is a positive integer indicating the mesh size and the level of polygonal approximation to the 
fractal boundary. With this choice of h, there are -/Vnlr(^) = 1 + (4 • 9 — 9 ■ 4 )/5 grid points in 
as well as 3 • 4^ grid points on dQ (see Table [1]). The zero-Dirichlet boundary conditions are 
imposed by setting m = at the grid points on the boundary. 

We used a different triangular grid. We found more accurate results with a larger h by choosing 
the grid spacing to be h = ^nss(^) =2/3 and placing the boundary between grid points. This 
yields N = Nnss(£) = ~~ 4 )/5 grid points in the snowflake region Q. No grid points are on 
<9f2 with our choice. We use ghost points, which are outside the region, to enforce the boundary 
conditions, as described in Section [2l To compare our results with those of ILapidus et al., 1996] 

u) th W) 

we will use A^ to denote the k eigenvalue of L at level £ with our method, and [x k to denote 
the eigenvalues published in |Lapidus et al., 19 96 . 



£ 


1 


2 


3 


4 


5 


6 


N NSS (£) 


1 


13 


133 


1261 


11605 


105469 


N LNR (£) 


1 


37 


469 


4789 


45397 


417781 



Table 1. The number of interior grid points with h = /inss(^) = 2/3 (the current 
work), and with h = /ilnrCO = 1/3^ as in Lapidu s et al., 1996| . Note that our 
new grid has approximately 75% fewer grid points at the same £, when £ is large. 
The region fi is open and the larger values of TV" published in |Lapidus et al., 1996 
includes the grid points on the boundary. 

A different approach, avoiding triangular grids altogether, can be found in the unpublished thesis 
|Banjai, 2003| . This work uses the conformal mappings found in |Banjai & Trefethen, 2003 . 

In |Neuberger and Swift, 200"T| , the Gradient Newton Galerkin Algorithm (GNGA) was devel- 
oped to investigate existence, multiplicity, nodal structure, bifurcation, and symmetry of problems 
of the form (I5.1D . (This PDE is found in the concluding section.) The GNGA requires as input 
an orthonormal basis of a sufficiently large subspace consisting of eigenfunctions of the Laplacian. 
Since our eventual application concerns solving the nonlinear equation (I5.1h on the region Q with 
fractal boundary, we face the considerable challenge of first obtaining eigenfunctions (solutions to 
the linear problem (11. ip ) numerically. In |Lapidus et al., 1996| , this was done using essentially the 
inverse power method with deflation on a triangular grid. They used approximating boundary 
polygons with vertices at grid points. 

Using our new grid, we improve upon their results and are substantially successful in obtaining 
a basis of such functions for a sufficiently large subspace for our future nonlinear needs. We use 
the sophisticated numerical package ARPACK instead of deflation. This software is based upon 
an algorithmic variant of the Arnoldi process called the Implicitly Restarted Arnoldi Method (see 
|Lehoucq et al., 1998| ) and is ideally suited for finding the eigen-pairs of the large sparse matrices 
associated with the discretization L of the Laplacian. It is easily implemented, requiring only a 
user-provided subroutine giving the action of the linear map. One of our innovations in investigating 
the snowflake is taking the boundary to lie between grid points and using ghost points just outside 
f2 when approximating the Laplacian at interior points closest to the boundary. This results 
in better approximations of true eigenvalues using fewer interior grid points than achieved by 
|Lapidus et al., 19 96] using the standard grid method of enforcing the boundary condition. We 
support this claim by comparing our results via curve fitting data points to predict the true values. 

In Section 2, we describe in more detail the triangular grid and the accompanying second dif- 
ference scheme for approximating the Laplacian, as well as the ARPACK implementation using 
this information to generate the basis of eigenfunctions. In Section 3, we compare our numerical 
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eigenvalue approximations to those obtained in Lapid us et al., 1996| . In particular, we perform 
Richardson extrapolations on both data sets. In Section 4, we apply representation theory to de- 
termine the 8 possible symmetries that eigenfunctions (and approximating eigenvectors) can have, 
given the H)q symmetry of the region f2 and the approximating grids. We consider this rigorous 
treatment of symmetry to be a key contribution of this paper. This information is used for numerical 
post-processing to find symmetric "canonical" representatives for multiple eigenvalues. This cata- 
log of the symmetries of basis elements will be used in an essential way in our subsequent nonlinear 
bifurcation studies. Section 4 also contains graphics depicting a selection of approximating eigen- 
vectors to both problems in (jl.ip . Section 5 gives a brief indication of how our new grid can be used 
when implementing GNGA on related nonlinear problems (see [Ncubcrg er, Sieben, and Swift II| ). 
Also, we discuss how the known symmetries can be exploited to reduce the number of integrations 
required by that scheme. 

2. Ghost Points and ARPACK. 

In approximating the Laplacian for functions defined on 0, we developed the grid technique 
depicted in Figure [TJ As in [Lapidus et al., 1996| , at interior points with interior point neighbors 
one sees that the standard second-differencing scheme when applied to a triangular grid leads to 
the approximation 

—Au(x) « —-g ^Qu(x) — ^^{6 neighbor values of u}^ . 

Our scheme differs, however, when computing approximations at interior grid points with neighbors 
that lie outside the boundary. In |Lapidus et al., 1996| the value of zero at boundary points (which 
lie on their grid) is used to enforce the zero-Dirichlet boundary condition. When we approximate the 
Laplacian at a point X{ near the boundary, we set u = —u(xi) at ghost points which are neighbors 
of x\. Specifically, for the level I = 2 example found in Figure [JJ and with the understanding that 
Ui ~ u(xi), we have 
2 

-Au(xi) « (Qui - (u 2 + u 3 + 1i 4 + u 5 + u e + u 7 )), 

2 2 

-Au(x 2 ) « ^2"( 6n 2 - ((-U-2) + Ui+U 3 + U 7 + U 8 + tig)) = ^V u 2 ~ ( u l + u 3 + u 7 + u 8 + Ug)), 

2 2 

-Au(xg) « ^2"( 6n 9 - ((~ u 9) + (-Ug) + (-Ug) + (-Ug) + U 2 + tt 3 )) = ^ (lOtig - (u 2 + ti 3 )). 

In the first line, there are no ghost points used because all the neighbors of xi are interior points. In 
the approximation at x 2 , (—u 2 ) represents the value of u at the ghost point gi, as labelled in Figure 
[TJ In the last line, (—tig) represents the value of ti at gi, . . . ,54. Note that the value of u at gi is 
different in the calculation at x 2 and xg. An alternative way of imposing the boundary conditions is 
to set ti((?i) to be the average of u 2 , us, and ug. We did experiments with this alternative method, 
but the results were not as accurate as the method we have described. 

Our method of imposing the zero-Dirichlet boundary conditions can be summarized as 

— Au(x) ~ — —^((12 — (number of interior neighbors))u(a;) — ^^{interior neighbor values of u}). 

The zero-Neumann boundary condition can be easily applied as well. Indeed, one has the only 
slightly different formula (which agrees away from the boundary) 

2 x ^ 

—Au(x) ~ ^-^((number of interior neighbors) u(x) — ^^{interior neighbor values of u}). 

To enforce the Neumann instead of the Dirichlet condition, we only need to change 3 characters of 
our ARPACK code. Specifically, one deletes the "12—" found in the zero-Dirichlet formulae above. 
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Figure 1. The Koch snowflake <9f2 with iV N ss(2) = 13 labelled grid points {xi}^ 
at level t = 2. The grid used by |Lapidus et al., 1996 consists of the A ? lnr(2) = 37 
large and small points inside the snowflake, along with 48 small points on <9f2. The 
points outside of the snowflake, some labelled g^ are ghost points we use to enforce 
the boundary conditions. For example, u(g-i) = —u(xg) for Dirichlet boundary 
conditions and u(g^) = u(xg) for Neumann boundary conditions. On the other 
hand, u takes on different values at g\ when the Laplacian is evaluated at x 2 > x 8> or 
xg. 



The user-provided ARPACK subroutine takes as input a vector v G M. N and outputs w G W N 
with w = Lv for the N x N matrix L approximating the discretized negative Laplacian, where 
for convenience we use N to denote N^ssi?)- This procedure is easily coded once an JV x 6 
dimensional array t with neighbor information is populated. In the pseudocode in Figure[2l t(i,j) G 
{0, 1, 2, ... , N} is the index of the j th neighbor of the grid point x%, j = 1, . . . , 6. If t(i,j) = for 
some j, then the i th interior point is near the boundary and has less than 6 interior point neighbors. 
We let ki G {2, 4, 5, 6} denote the number of interior point neighbors of grid point X{. 



wit) 



for Dirichlet boundary conditions, or 
for Neumann 



Loop for i = 1, . . . , N 
1. Set 

(12 - ki) * 

ki * v(i) 

Loop for j = 1, . . . , 6 

a. Find index p = t(i,j) of j th neighbor 

b. If p ^ then subtract neighbor value: w(i 
Multiply by h factor: w(i) = 2 * w(i)/ (3.0 * h * h) 



wit) 



v(p) 



Figure 2. Pseudo code for user-provided subroutine encoding the linear map v i— > 
w = Lv. 



The neighbor file is generated using the set and vector data structures and the binary search 
algorithm of the Standard Template Library in C++. In the first step we find the integer coordinates 
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of the grid points in the basis {(1,0), 2 )}• ^ ne P roce dure uses simple loops to find the coor- 
dinates of the grid points inside a large triangle and then calls itself recursively on three smaller 
triangles placed on the three sides of the original triangle, until the desired level is reached. To 
avoid duplication of grid points the coordinates are collected in a set data structure. In the second 
step, we copy the coordinates into a vector data structure and use binary searches to find the 
indices of the six possible neighbors of each grid point. In the last step, we compute the Cartesian 
coordinates of the grid points and write them into a file together with the indices of the neighbors. 

3. Numerical Results 

In this section we present our experimental results. Our best approximations A^ for the eigen- 
values are obtained by performing Richardson extrapolation. Specifically, we find the y-intercepts 
of the Lagrange polynomials fitting the points {(/inssCO) ^*(^))}|=4- We also compute Richardson 
extrapolations fj^ using the data published in |Lapidus et al., 1996| . In Table [2] we list the level 6 
and Richardson approximations of the first ten and the 100 th eigenvalues. 





NSS 


NSS 


LNR 


LNR 


k 


A fc (6) 




M*(6) 


$ 


1 


39.353 


39.349 


39.390 


39.352 


2 


97.446 


97.438 


97.537 


97.438 


3 


97.446 


97.438 


97.537 


97.438 


4 


165.417 


165.409 


165.622 


165.478 


5 


165.417 


165.409 


165.622 


165.478 


6 


190.381 


190.373 


190.571 


190.365 


7 


208.622 


208.617 


208.837 


208.59 


8 


272.415 


272.413 


272.755 


272.480 


9 


272.415 


272.413 


272.755 


272.480 


10 


312.348 


312.358 


312.645 


312.351 


100 


2322.129 


2324.925 







Table 2. The first ten and the 100 eigenvalues to the Dirichlet problem. NSS 
denotes our new grid, while LNR denotes the grid in [Lapidus et al., 1996| . Provided 
are the level £ = 6 approximations, where the NSS scheme uses A^nss(6) = 105469 
grid points and the LNR method uses ./Vlnr(6) = 417781 interior grid points. The 
symbols A? and //? denote the the Richardson extrapolation values for A& using 
levels I G {4, 5, 6}. Blank entries correspond to no data available for comparison. 



We computed the relative differences (Afe(6) — A^) /X^ and U£fc(6) — fi^j //x]^ for k G {1, . . . , 10}. 
The relative differences of the A values ranged from 10~ 5 to 10 -4 , while the relative differences 
of the values are all on the order of 10 -3 . The absolute differences between the Richardson 
extrapolations A^ and fj^ ranged from 10~ 4 to 7 • 10~ 2 . Note that even for k = 100 we have 
(Aioo(6) - X%) /X% « -1.203 • 10" 3 . 

In Figures [3] and U] we visually compare the Richardson extrapolations for Ai and Aio- We can 
see that although the extrapolated values are nearly identical, our approximations are much closer 
to the common extrapolated values using a lot fewer grid points. This is a key issue for us since 
in |Neuberger, Sieben, and Swift II| we will require accurate eigenvectors and eigenvalues using as 
few grid points as possible. We cannot use the extrapolated eigenvalues, since although they are 
more accurate approximations of eigenvalues of —A, they are not eigenvalues of L corresponding 
to eigenvectors of L at any given level. 
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Figure 3. Two views of the Richardson extrapolations for Ai. The solid line is the 
graph of the Lagrange polynomial fitting our data for I e {3,4,5,6}. The dashed 
line fits the i G {3, 4, 5, 6} data of |Lapidus et al., 1996| together with the unpub- 
lished level i = 7 eigenvalue approximation obtained via private communication 
from Robert Renka. 
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2/3 4 



Figure 4. Two views of the Richardson Extrapolations for Aio- As in Figure El 
the solid lines correspond to our data and the dashed lines to the data in 
|Lapidus et al., 1996| . 



Our results for the Neumann boundary conditions are shown in Table [3l Based on private com- 
munication we know that our approximate eigenvalues are very close to the unpublished numbers 
obtained by Lapidus et al. A careful comparison of the two grid schemes is not possible at this 
time, since we do not have all of their data. We found that the Lagrange polynomial of our Neu- 
mann data approaches h = linearly, similar to the curves of |Lapidus et al., 1996| in Figures [3] 
and|H This led us to consider an alternate scheme for enforcing boundary conditions. As noted in 
Figure [Q u is multi- valued at certain ghost points. We tried using a single average value at these 
ghost points, but the slope of the Lagrange polynomial at was larger, so that the eigenvalues at a 
given level were farther from the extrapolated value. It is an area for future research to understand 
why the ghost points work so well. All we can now assert is that our method clearly out-performs 
that found in |Lapidus et al., 1996| for the zero-Dirichlet eigenvalue problem on the snowflake re- 
gion. The ghost points can be used in general regions, and it would be interesting to determine the 
optimal method for enforcing boundary conditions on general regions. 

We produced contour plots of the eigenfunctions. An example is shown in Figure [5j These 
contour plots were produced by a Mathematica notebook that reads in the u vector and outputs 
a postscript file. The level of the grid approximation is computed from the length of the u vector. 
All of the contour plots shown in this paper use I = 5 data for which the u vector has length 
Wnss(5) = 11605. 
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k 


A fc (6) 




(Afc(6)-A*)/A* 


1 


0.0000 


0.0000 


NA 


2 


11.9105 


11.8424 


0.0057 


3 


11.9105 


11.8424 


0.0057 


4 


23.1770 


23.0466 


0.0057 


5 


23.1770 


23.0466 


0.0057 


6 


27.5770 


27.4261 


0.0055 


7 


52.4164 


52.2105 


0.0039 


8 


85.8449 


85.5521 


0.0034 


9 


85.8449 


85.5521 


0.0034 


10 


112.7801 


112.0200 


0.0068 


100 


1295.4431 


1271.1900 


0.0191 



Table 3. The first 10 and 100* eigenvalues for the Neumann problem. We have 
included the level £ = 6 approximations, the Richardson extrapolations using £ G 
{4,5,6}, and the relative differences of the two. All results use our new grid. The 
Neumann eigenvalues obtained by Lapidus et al. have not been published. 




Figure 5. The graph, and a contour plot, of the sixth eigenfunction, ipg, with 
Dirichlet boundary conditions. The graph uses 469 grid points to triangulate the 
snowflake region. The contour plot shows our level £ = 5 data with the new grid, 
as described in the text. The white and black regions in the contour plot represent 
positive and negative values of ipQ. The contours are equally spaced, and the dots 
represent local extrema of ipQ. 

The local extrema of u are calculated in two steps. First the extreme values of U{ are calculated. 
Then, a quadratic fit to this data point and its six neighbors is performed. A dot is then drawn 
at the extremum of the quadratic function. This extra effort, compared to drawing a dot at the 
grid point, has a noticeable effect even at level 5. After the extrema are found, the u values of the 
contours are computed using a heuristic that gives fewer levels as the number of extrema increases. 

The black regions are then drawn by subdividing the snowflake region into the triangles defined 
by the grid points. If u < at all three vertices of a triangle then the triangle is filled with black. 
If u < at some vertices of the triangle and u > at others, then a linear interpolation is used to 
estimate the region where u < within the triangle. The contours are also produced using linear 
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interpolation within the triangles: If the value of u on the vertices spans a contour value, then a 
short line segment inside the triangle is drawn based on the linear fit. 

Several details of the implementation of the contour plotting have been left out. For example 
the region outside of the triangulation of grid points, but inside the snowflake boundary, is shaded 
by a different technique. 



4. Symmetry and the Canonical Basis. 

Some of the eigenvalues of the Laplacian on the snowflake, (jl.ip . have multiplicity one, and some 
have multiplicity two. In this section we quote well-known results in group representation theory 
to explain the observed multiplicity. We also describe a canonical way to choose two eigenvectors 
to span the two-dimensional eigenspaces. Details of group representation theory can be found in 
Tinkham, 1964] , [Sternberg, 19941 and |Scott, 1964| 



Assume that G is a finite group. A linear representation of G is a homomorphism a : G — > GL(U) 
where GL(U) is the group of invertible linear operators on the vector space U = R N or C^. The 
vector space U is called the representation space of the linear representation. If B is a basis for 
U and T G GL{U) then we write [T]b for the matrix of T in the basis B or simply [T] if B is 
the standard basis. We call the map g i— > \a(g)\ a matrix representation. Two representations 
a, j3 : G — )• H are equivalent, and we write a ~ (3, if there is an h G H such that j3(g) = h~ 1 a(g)h 
for all g G G. 

Let a : G — > GL(U) be a linear representation. If g G G then a(g) : U — > U is a linear operator; 
for convenience we sometimes use the notation a g = a(g). The linear representation a induces 
a group action G x U —> U . We often write g ■ u in place of oc g [u) when the representation a is 
understood. 

A subspace W of U is called an invariant subspace of a if a g (W) C W for all g £ G. The 
representation a is called irreducible if a has no proper invariant subspaces. The property of 
complete reducibility, also known as Maschke's theorem, says that there are a-invariant subspaces 
Ui, . . . , Uk such that U = U\ © • • • © Uk and 7 <n ) := a\jj n is irreducible for each n G {1, . . . , k}. If 
B n is a basis for U n and -B = L) n B n then the matrix of a g in the basis B is block diagonal for all 
g G G, that is 

fc 

[a 9 b = 0[7 (n) G?)b„. 

n=l 

Let rW, i G {1, be an element from each of the q equivalence classes of irreducible 

representations of G. Suppose we have a complete decomposition of the representation a into 
irreducible representations 7'"'. For each i there is an a-invariant subspace 

V d) = 0{C/ n I 7 W „ r (0}. 

Whereas there is great freedom in choosing the elements U n in U = fTiffi- • -ffif/jt, the decomposition 
[/ = yW © ... © y(?) 

is unique up to ordering. 
The characters of the representation are X^Hd) = tr[T®(p)]. These characters are used in 
projection operators 

(4.1) 

onto the invariant subspaces V® = (U) where di is the dimension of the i th irreducible repre- 
sentation. Note that [rW(g)] is a d{ x di matrix. 

To proceed we must choose a fixed set of matrix representations [T^] from each equivalence 
class of irreducible representations. We call these canonical matrices. The following calculations 
are simplified if we make the canonical matrices as simple as possible. It is always possible to 
choose matrices that are unitary, and we assume that the canonical matrices are unitary. 
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There is set of projection operators described in |Tinkham, 1964| 
(4-2) P^ = ^Ei T(t) (9)k 3 » g , 

where the coefficient of a g is the j th diagonal element of the matrix [r^ (<?)]. Note that = 

Sfci^j • The corresponding vector spaces vj = Pj l \u) are not a-invariant, but they are 
orthogonal, and we shall see that the following decomposition of the representation space is very 
useful: 

»=i \j=i 

We now consider the effect of the symmetry on commuting linear operators. The theory is much 
simpler if we assume that the representation space U is complex. 

Schur's Lemma. Suppose T\ and T2 are two irreducible representations of G on C dl and C d2 , 
respectively, and S : C^ 1 —¥ C d2 is a linear operator such that ST\(g) = T2(g)S for all g G G. Then 
S = ifV\ and T2 are not equivalent, and S = cl for some c € C ifY\ = T2. 

Note that Schur's Lemma does not address the case where Ti and T2 are different, but equivalent, 
representations. This case can be addressed by choosing a basis where Ti and T2 are the same, 
as in the proof of the following corollary. An abstract version of this corollary can be found in 
|Sagan, 1991] Theorem 1.2.8]. 

Corollary 4.4. Suppose that U = C N is a representation space for a that is decomposed as in 
$i4-3ty , and T : U —> U is a linear operator that commutes with a. That is, a g T = Ta g for all 

g £ G. Then each of the spaces vj is T -invariant, and the operators Tj := T\ (*) decompose T 
as 

(4.5) T = ©(©!T< 

i=i \j=i 

Furthermore, Tj is similar to Tj, . . 

Proof. By complete reducibility, U = © n= i U n , where a\u n is equivalent to the irreducible repre- 
sentation r^, with i n G {l,...,g}. We can write T in terms of k 2 blocks T m , n — P\j m o T\y n : 
U n — > U m . For each n, we can choose a basis U n = span{e nj - | j = l,...,di n } such that 

Pj e n,j' = 3i,i n dj,j' e n,j, and a g (e n j) = Yl^=i e n,j'[^ in \9)]j' ,j ■ I n other words, the basis vectors 
e n j G U n transform like the standard basis vectors ej G IR di ™ do when multiplied by the matrices 
[r( l ")(g)]. If i n = i m , then a\u n = a\u m - Thus, Schur's Lemma implies 



m=l 

where c m ^ n = if i m ^ i n . So, if i m = i n then T(e n j) is a sum over e m j with the same j. The 

of projection operators ca 

vj 1 ' = span{e nj - | i n = i}, 

lvariant. Finally, Tj is si 
basis we have constructed they are represented by the same matrix. □ 



m=l 

= if i m ^ i n . So, if i m = i n thei 

(i) 

spaces Vj that we defined in terms of projection operators can be written as 

r(t) 
j 

and it is now clear that vj^ is T- invariant. Finally, is similar to Tjj ^ if i = i', since in the 
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Remark 4.6. The spectrum of T is the set of eigenvalues of the k x k matrix C, with elements 
Cm,n- The matrix C separates into q diagonal blocks one for each irreducible representation 
class of G. If A is an eigenvalue of then A is an eigenvalue of T with multiplicity at least 
di- If the matrix C has any multiple eigenvalues, then we say that T has accidental degeneracy. 
In this case, a perturbation of T that commutes with a can be chosen so that the perturbed C 
matrix has simple eigenvalues. Barring accidental degeneracy, every eigenvalue of T corresponds 
to a unique irreducible representation r'*), the eigenvalue has multiplicity di, and an orthonormal 
set of eigenvectors of this eigenvalue can be chosen, one from each of the di subspaces Vj , j £ 
{1,2,..., di}. These d, eigenvectors can be used as a basis for one irreducible space U n . In this 
way, we can write U = ® n= i U n , where each irreducible component U n is the eigenspace of an 
eigenvalue A n of T. 

If we require the spectrum of a linear operator T : W N — > M. N that commutes with a representation 
on a real vector space, we first complexify to the operator T : C N — > <C N . The spectra of T and 
T are the same, so the consequences of Schur's Lemma hold provided we consider representations 
rW which are irreducible over the field C. Sometimes an irreducible representation over R breaks 
up into two complex conjugate irreducible representations over C. The eigenvalues of T associated 
with this pair of representations are complex conjugates. If T is self-adjoint, and thus has real 
eigenvalues, this causes "accidental degeneracies." For the Laplacian operator on the snowflake 
domain we do not have to consider this complication since the irreducible representations we need 
are the same over the field of real or complex numbers. 

We now apply the general theory to the eigenvalue equation (jl.lj) . The symmetry group of the 
snowflake region, as well as the set of grid points in f2, is the dihedral group 

D 6 = (p, a | p 6 = a 2 = I, pa = ap 5 ). 

It is convenient to define r = p 3 a. Note that ar = ra = p 3 , and p 3 commutes with every element 
in Dq. The standard action of IS)q on the plane is 

p ■ fa v) = (h x + ~^t x + b) 

(4-7) a-(x,y) = (-x,y) 

t ■ (x,y) = {x,-y). 

In this action p is a rotation by 60°, a is a reflection across the y-axis, and r is a reflection across 
the x-axis. 

For a given grid with iV points in fi, the ©6 action on the plane (|4,7p induces a group action on 
the integers {1,2,3,..., N} defined by x g .% = g ■ Xj. There is also a natural action on the space of 
all functions from U to K, U = M. N , defined by (g ■ u){xi) = u(g _1 ■ Xi) for all u G U and g £ ©6- 
With the usual identification it, = u(xi), this action corresponds to a linear representation a of H)q 
on U = 1^ defined by (a g (u))i = u g -i.i. 

There are exactly 6 irreducible representations of Dg up to equivalence. Our canonical matrices 
are listed in Table HI Since we have chosen a set of real matrices, and P- are projection 
operators on both M. N and C^. We can write the representation space U = M. N as 

u = v« e e e e v} 5) e e v} 6) e y 2 (6) , 

where each of the real vector spaces V W and Vj is invariant under any commuting linear operator 
on U. There are just six spaces in this decomposition regardless of the number of grid points N. 
On the other hand, the decomposition U = J2n=i has thousands of irreducible components U n 
when iV is large. 

In the remainder of this section we use T to refer to the negative Laplacian operator. The 
results of Corollary 14.41 hold, since T commutes with the action on U. We will use the terms 
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i 


[rW( P )] 


[rW(a)] 


[rW(r)] 


1 


i 


1 


1 


2 


i 


-1 


-1 


3 


-l 


1 


-1 


4 


-l 






K 


/ -1/2 V3/2 \ 

V -V5/2 -1/2 y< 


/I \ 


/I \ 
[o -l) 


6 


/ 1/2 \ 
V -V3/2 1/2 y 


(o -l) 


(o?) 



Table 4. A representative from each of the six equivalence classes of irreducible 
matrix representations of J}q. These canonical matrices are real, but the represen- 
tations are nevertheless irreducible over the complex numbers. The homomorphism 
condition T^\gh) = T^\g)T^ (h) allows all of the matrices to be computed from 
[rw(p)] and [rW(p)]. The last column is included for convenience. The representa- 
tion r^ 6 ) corresponds to the standard action of H>6 on the plane (14. 7h . 

eigenvectors and eigenfunctions of T interchangeably, since the vector u G is a function on the 
iV grid points. 

The numbering, 1 through 6, of the irreducible representations of Dq is somewhat arbitrary. 

(i) 

Therefore we give the names V PxPy( i to the 8 T-invariant spaces V- , where p x and p y describe the 
parity of the functions under the reflections a and r respectively, and d is the dimension of the 
associated irreducible representation. Some calculations allow us to give simple descriptions these 
T-invariant spaces without having to appeal to the projection operators: 



V ++ i 


:= 


= {u£U | p 


u = u, a ■ u = u, t ■ u = u} 




V—x 


:= V {2) 


= {u e U | p 


u = u, a ■ u= —u, t • u = - 


-u} 


V+-i 


:= 


= {u G U | p 


u = —u, a ■ u = u, t • u = - 


-u} 


V-+i 


:= 


= {u£U | p 


u = —u, a ■ u = —u, t ■ u = 


«} 






= {u£U\p 3 


■ u = u, u + p 2 ■ u + p 4 ■ u = 


0} 






= {u€U\p 3 


■ u = —u, u + p 2 ■ u + p 4 ■ u 


= 0} 


V++2 


:= V} 5) 


= {u G V {5) | 


a ■ u = u, t ■ u = u} 




V—2 


— 

■— v 2 


= {u g y (5) | 


a ■ u = —u, t ■ u = —u} 




V+-2 


:= V} 6) 


= {u g y (6) | 


a ■ u = u, t ■ u = —u} 




V-+2 


— 

■— v 2 


= {u g y (6) | 


a ■ u = —u, t ■ u = u} 





We can use Corollary 14.41 in two different ways. First, we could construct the numerical approxi- 

mations to the restricted operators T- and find the eigenvalues and eigenvectors of the correspond- 
ed) 

ing matrices. This has the advantage that the matrices which represent are smaller than those 
representing T. In section [2] we described the simple structure of the matrices L, which represent 
T in the standard basis. We have not attempted to construct the more complicated matrices which 
represent T- . 

We use a second approach. We compute eigenvalues and eigenvectors of the full operator T, 
and use Corollary 14.41 to classify the eigenvalues according to the corresponding irreducible rep- 
resentation. An eigenvector tp of T will satisfy P'''(^) = ip for exactly one i, barring accidental 
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degeneracy. By computing all the projections, we can determine which irreducible representation is 
associated with the eigenvector. If the irreducible representation has dimension greater than one, 

then we can choose an orthonormal set of eigenvectors such that each eigenvector lies in a unique 

(i) (i) 
Vj using the projections P- . 

For our group G = 3q, we do not have to use the full projection operators pW and Pj % \ We 
only need the four projection operators 

_ qi + p x a a + p y a T + p x p y a p 3 
^PxPy — 7 > 

where p x ,p y £ {+, — }• For example P^ (u) = (u + a ■ u — r • u — /? 3 ■ u)/A. Note that P PxPy (U) = 

Vp x pyi © Vp x p y 2- If ip is an eigenvector with multiplicity 1 and P PxPy ip ^ 0, then P PxPy ip = ip and 
ip £ VpxPyl- ^ 1S noteworthy that we do not have to check for rotational symmetry. Similarly, if ip 
is an eigenvector with multiplicity 2 and P PxPy ip / 0, then P PxPy ip £ V PxP 2- 

Here is our algorithm to compute the symmetry of the eigenfunctions of T, and to find eigen- 

(i) 

vectors in the T-invariant spaces Vf'. First of all, we know from the computed eigenvalues which 
are simple and which are double. If an eigenvalue A is simple, then the corresponding eigenvector 
ip satisfies P PxPy ip = tp for exactly one choice of p x and p y , and ip £ V PxP \. If an eigenvalue A 
has multiplicity 2, then each of the 2 orthogonal eigenfunctions returned by ARPACK is replaced 
with the largest of the four projections P PxPy ip- Then, the projected eigenfunctions are normal- 
ized. If necessary, the eigenfunctions are reordered so that ip £ V++2 comes before ip G V 2, and 

tp £ V-\ 2 comes before ip E V-+2- This algorithm assumes that there is no accidental degeneracy. 

We produced a list of the eigenvalues and symmetry types up to the 300 th eigenvalue, confirming 
that there is no accidental degeneracy. 

We now explain why replacing each eigenfunction with the largest projection gives two linearly 
independent eigenfunctions. Suppose ip and <p are two eigenfunctions with the same eigenvalue. 

Then either ip, <p € = V++2 © V 2 or ip, (f> G = 2 © V-+2- Assume that we are in the 

first case. The eigenfunctions are orthonormal, so that H^H 2 = H^H 2 = 1 and the inner product 

(ip,4>) = 0. We want to replace ip and (j) with normalized eigenfunctions ip++ £ V++2 and ip £ 

V — 2- We can write ip = a^,ip ++ + b^ip — and <p> = a<pip ++ + b^ip — , where ||V>|| 2 = + = 1, 
\\<p\\ 2 = a| + 6^ = 1, and (ip,<p) = av a </> + b^bj, = 0. In other words, (a^,6^) and (a^,b^) are two 
orthogonal unit vectors in the a-b plane. Therefore, a\ > 65, if and only if 6^ < a?. Neglecting the 
numerical coincidence of a\ = tfy, this implies that replacing ip and cp by their largest projection will 
include eigenfunctions of both symmetry types. A similar argument holds for pairs of eigenfunctions 
in V + - 2 © V-+2- 

Figures [6] and [7] shows eigenfunctions with each of the 8 symmetry types for the Dirichlet and 
Neumann eigenvalue problems, respectively. More eigenfunction contour plots can be found at 
http : //jan.ucc .nau. edu/~ns46/snow/. 

5. Conclusions. 

Table [5] summarizes our best estimates of the eigenvalues, and the symmetry types of the eigen- 
functions, for the first few eigenvalues of the linear problem (jl.ip . We have no error bounds but 
we can get an indication of the accuracy by comparing our results with those in Banjai's thesis 
Banjai, 2003;. Banjai computed the first 20 eigenvalues of the Dirichlet problem, and our eigen- 
values agree with his results to four significant figures. This suggests that our results are accurate 
to four significant figures. Banjai's results for the first ten Dirichlet eigenvalues are more precise. 

We are ultimately interested in the connections between the linear problem (jl.ip and superlinear 
elliptic boundary value problems of the form 



(5.1) 



Au + f(u) = in n 

u = on (90. 
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Figure 6. The first occurrences of the 8 symmetry types of eigenfunctions of the 
Laplacian with zero Dirichlet boundary conditions. Each of the eigenfunctions in 
the first row has a simple eigenvalue. Those in the second row come in pairs. 
For example, A4 = A5, since ip4 and "05 are each in the invariant subspace 
corresponding to the 2-dimensional irreducible representation . Similarly, A2 = 
A3, since ip2 and tp3 are each in The eigenfunctions shown respect the canonical 
decompositions = v} 5) K, (5) and = v{ 6) V 2 (6) . 

The so-called Gradient-Newton-Galerkin- Algorithm (GNGA, see [NS]) seeks approximate solutions 
u = YljLi a ji )t j to (|5.ip by applying Newton's method to the eigenfunction expansion coefficients 
of the gradient VJ(it) of a nonlinear functional J whose critical points are the desired solutions. 
In [NSS2], we will enumerate the 23 isotropy subgroups of H> e x Z 2 , along with the associated fixed 
point subspaces. These are the possible symmetry types of solutions to (15. ip when / is odd. As 
a result, we will be able to follow the bifurcation branches of (|5.ip by varying the parameter A in 
the nonlinearity defined by f(u) = Xu + u 3 . Our symmetry information is crucial to the branch 
continuation decision making process. In [NSS2] we will find solutions of all 23 symmetry types. 

Our catalog of symmetry information is helpful in another way. Given that our basis is chosen in 
a symmetric fashion and by knowing which invariant subspace a given solution branch's elements 
lie, it is often possible to know that many of the eigenfunction expansion coefficients are zero. 
The GNGA requires numerical integration using values at the N = Anss(^) grid points for each 
coefficient of the gradient VJ(u), as well as for each entry of the Hessian matrix representing 
D-zJ{u). If we use a basis of eigenfunctions spanning a subspace of dimension M, this amounts to 
M 2 + M integrations. Using symmetry, we are able to avoid many of these integrations, resulting in 
a substantial speedup. This becomes very important as we seek high-energy solutions with complex 
nodal structure, where large values of M and N are required. Each solution is represented by a 
point on a bifurcation curve, whereby many solutions are sought to complete the branch. Many 
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Figure 7. The second occurrences of the 8 symmetry types of eigenfunctions of the 
Laplacian, with Neumann boundary conditions. Compare with Figure [6l 



of the branches of high-energy solutions have a rich bifurcation structure, resulting in multiple 
secondary and tertiary branches. There is a substantial time savings obtained by cutting down on 
the number of required integrations for each Newton step for each solution, on each branch. 

In conclusion, this paper presents an efficient, accurate, and easy to implement method for 
obtaining a basis of eigenfunctions to a subspace which is sufficiently large for performing the 
eigenfunction expansions required by the GNGA method in solving nonlinear problems of the form 
()5.ip . The major innovations of this work are in the new way we enforce boundary conditions 
and in the way we are able to choose symmetric representatives of eigenfunctions corresponding 
to multiple eigenvalues. Our new techniques for generating contour plots are also noteworthy. By 
further decomposing function space according to symmetry, our nonlinear experiments are much 
more successful (see |Neuberger, Sieben, and Swift II| ). We are currently porting the code to a small 
cluster using the parallel implementation PARPACK. In the future, we expect to use the techniques 
of this paper to generate larger numbers of more accurate eigenvalue/eigenfunction approximations, 
and in so doing, be able to investigate more complicated phenomena and/or problems with regions 
Q in M. n , for dimensions n = 3 and higher. 
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k 


(°) 


Symmetry 


Ai 1 (N) 


Symmetry 


1 
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+ + 1 





+ + 1 


2 
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+ - 2 
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+ - 2 


3 
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- + 2 


11.842 


- + 2 


4 


165.41 


+ + 2 


23.047 


+ + 2 


5 


165.41 


- - 2 


23.047 


- - 2 


6 


190.37 
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+ - 1 


7 


208.62 


+ - 1 


52.210 
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8 


272.41 


+ - 2 


85.552 


+ - 2 


9 


272.41 
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10 


312.36 
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112.02 
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11 


314.45 
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112.02 
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12 
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2 
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13 


359.53 
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14 


425.40 
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15 


443.54 
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139.54 
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16 
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- + 2 
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17 


458.66 
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151.00 
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18 
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19 


560.41 
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20 


560.41 


- + 2 


197.50 
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21 


566.79 


+ + 2 


207.10 


+ -2 


22 


566.79 


- - 2 


207.10 


- + 2 


23 


595.18 


+ + 1 


216.81 


+ + 2 


24 


617.70 


- - 1 


216.81 


- - 2 



Table 5. The first 24 eigenvalues of the negative Laplacian, and the symmetry of 
the corresponding eigenfunctions, with Dirichlet or Neumann boundary conditions. 
The Richardson extrapolations for the eigenvalues are given. A comparison with 
Banjai's thesis suggests that our results are accurate to four significant figures. The 
symmetry indicates the space containing the eigenfunction. For example, the first 
eigenfunction is in V++i for both boundary conditions. We follow the convention 
that V ++ 2 is listed before V. 2, and 2 is listed before V-+2- 



[Lapidus et al., 1996] M .L. Lapidus, J. W. Neuberger, R. L. Renka, and C. A. Griffith, Snowflake Harmonics and 

Computer Graphics: Numerical Computation of Spectra on Fractal Drums, International Journal Bifurcation and 

Chaos 6 no. 7 (1996), pp 1185-1210. 
[Lapidus and Pang, 1995] M. L. Lapidus and M. M. H. Pang, Eigenfunctions on the Koch Snowflake Domain, Com- 

mun. Math. Physics 172 no. 2 (1995), pp 359-376. 
[Johnson and Riess, 1982] L. Johnson and R. Riess, Numerical Analysis, Addison- Wesley : Reading, Mass. (1982). 
[Lehoucq et al., 1998] R. B. Lehoucq, D. C. Sorensen, and C. Yang ARPACK users' guide: Solution of large-scale 

eigenvalue problems with implicitly restarted Arnoldi methods. Software, Environments, and Tools. Society for 

Industrial and Applied Mathematics (SIAM) : Philadelphia, PA. (1998). 
[Neuberger and Swift, 2001] John M. Neuberger and James W. Swift Newton's method and Morse index for semilinear 

elliptic PDEs International Journal Bifurcation and Chaos 11 no. 3 (2001), pp 801-820. 
[Neuberger, Sieben, and Swift II] John M. Neuberger, Nandor Sieben, and James W. Swift Symmetry of Solutions 

to Semilinear Dirichlet Problems on Koch's Snowflake, preprint, 2003. 
[Sagan, 1991] B. E. Sagan The symmetric group, Wadsworth & Brooks/Cole Adv. Books Software, Pacific Grove, 

CA, 1991. 

[Scott, 1964] W. R. Scott Group theory, Prentice Hall, Englewood Cliffs, N.J., 1964. 

[Sternberg, 1994] S. Sternberg, Group Theory and Physics Cambridge University Press, Cambridge (1994). 
[Tinkham, 1964] M. Tinkham, Group Theory and Quantum Mechanics McGraw-Hill Inc., New York (1964). 



16 JOHN M. NEUBERGER, NANDOR SIEBEN, AND JAMES W. SWIFT 

E-mail address: John.Neuberger@nau.edu, Nandor.Sieben@nau.edu, Jim.Swift@nau.edu 

Department of Mathematics and Statistics, Northern Arizona University PO Box 5717, Flagstaff, 
AZ 86011-5717, USA 



